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(57) ABSTRACT 

Each voxel is assumed to contain exactly two distinct 
materials, with the volume fraction of each material being 
iteratively calculated. According to the method, the spec- 
trum of the X-ray beam must be known, and the attenuation 
spectra of the materials in the object must be known, and be 
monotonically decreasing with increasing X-ray photon 
energy. Then, a volume fraction is estimated for the voxel, 
and the spectrum is iteratively calculated. 

8 Claims, 2 Drawing Sheets 
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METHOD FOR BEAM HARDENING 
CORRECTION IN QUANTITATIVE 

COMPUTED X-RAY TOMOGRAPHY 

RELATED APPLICATIONS 5 

The present application claims the benefit of priority from 
copending provisional patent application 60/108,257, filed 
on Nov. 12, 1998, which is hereby incorporated by refer- 
ence. 10 

This invention was made with U.S. Government support 
under grant numbers NCC-5186 and NCC2-5088 awarded 
by NASA. The National Institutes of Health also supported 
the development of the present invention under grant P41- 
RR09784. The Government has certain rights in the inven- ^ 
tion. 

FIELD OF THE INVENTION 

The present invention relates generally to computer 
assisted tomography and quantitative computed tomography 20 
(both known as CT imaging). More particularly, it relates to 
a method for processing CT imaging data that corrects for 
beam hardening errors. 

BACKGROUND OF THE INVENTION 25 

CT imaging is extensively used for medical imaging and 
the imaging of objects. In CT imaging, X-rays are projected 
through the object being imaged, and these X-rays are 
detected by arrays of detectors. The X-rays are projected 3Q 
through the object in many different directions. The com- 
bination of X-ray trajectories through the object provides 
data from which the internal structure of the object can be 
determined. Contrast in CT images is provided by variations 
in X-ray attenuation within the object. No other contrast 35 
parameters are available. Therefore, accurate measurements 
of X-ray attenuation are required for high quality CT 
images. 

Quantitative computed tomography (QCT) is a technique 
that allows for quantitative measurements of physical prop- 40 
erties related to X-ray attenuation. QCT has been used for 
in-vivo quantitative measurements of bone density, for 
example. The uses of QCT include assessment of spinal 
trabecular bone, evaluation of drug therapy in the treatment 
of osteoporosis, screening for osteoporosis, fracture risk 45 
assessment and many others. Although QCT is now an 
established tool for bone densitometry, there exist major 
issues affecting the accuracy and precision of QCT mea- 
surements. 

QCT measurements and CT images are affected by beam 50 
hardening error in X-ray attenuation measurements. Beam 
hardening error is caused by the energy-dependence of 
X-ray transmission within an object being imaged. In any 
material, low-energy X-rays are attenuated more strongly 
than high-energy X-rays. Therefore, as a polychromatic 55 
X-ray beam passes through an object, the proportion of high 
energy X-rays in the beam increases, and attenuation 
decreases. Long path lengths through an object therefore 
appear to have an excessively small attenuation. When an 
image is computed, the center of an object appears to have 60 
a lower attenuation than the outer regions of the object. In 
this way, beam hardening error produces inaccurate mea- 
surements in QCT. Correction of beam hardening errors has 
been an active area of research since 1975. Some popular 
current correction techniques for beam hardening require 65 
strict assumptions about the X-ray attenuation characteris- 
tics of the materials within the object. The two most com- 


2 

monly used techniques are the water and bone corrections 
which assume that the materials in the scan field are either 
water-equivalent or dense bone -equivalent in X-ray attenu- 
ation characteristics. For more information regarding these 
techniques, reference can be made to “A Method for Cor- 
recting Bone Induced Artifacts in Computed Tomography 
Scanners” by R M. Joseph and R. D. Spital in the Journal 
of Computer Assisted Tomography vol. 2, pp. 481-487, 1978 
and “Post-Reconstruction Method for Beam Hardening in 
Computed Tomography” by O. Nalcioglu and R. Y. Lou in 
Physics in Medicine and Biology, vol 24, pp. 330-340, 1979. 

In another beam hardening correction method, calibration 
tubes having a known transmission characteristic are used. 
However, this method is often not effective for two reasons: 

1) different regions in the scan field experience different 
degrees of beam hardening, and 2) calibration tubes cannot 
capture the beam hardening characteristics in vicinity of a 
patients bone because the calibration tube must be placed 
outside the body. 

There exists a need in the art of CT imaging and QCT for 
an improved method of beam hardening correction. An 
improved beam hardening correction technique will provide 
QCT measurements and CT images with improved accuracy. 

OBJECTS AND ADVANTAGES OF THE 
INVENTION 

Accordingly, it is a primary object of the present invention 
to provide an improved method for beam hardening correc- 
tion that: 

1) accurately corrects beam hardening errors; 

2) does not require knowledge of the attenuation character- 
istics of the X-ray detectors used in the CT imager; 

3) does not require calibration tubes; 

4) can be used with objects comprising many different 

materials. 

These and other objects and advantages will be apparent 
upon reading the following description and accompanying 
drawings. 

SUMMARY OF THE INVENTION 

The present invention provides a method for beam hard- 
ening correction in CT imaging data. The present method 
includes reiterative calculations that converge on accurate 
measurements of X-ray attenuation. In a first variation of the 
present invention, the following information is required: 

1) The attenuation spectra for materials within the object 
being imaged. Each voxel is assumed to contain at most 
two materials. 

2) The output spectra of the X-ray source. 

3) The output data from the X-ray detectors. 

An initial estimation is made of volume fraction of the 
two materials in each voxel. Then, a reiterative calculation 
is performed that converges upon the true volume fraction 
for each material in each voxel. 

In a second variation of the present invention, the fol- 
lowing information is required: 

1) Two basis attenuation spectra. The attenuation spectra 
do not need to correspond to real materials. 

2) Output spectra from the X-ray source at two different 
settings (e.g. two different X-ray tube voltages). 

3) Output data from the X-ray detectors at the two X-ray 
source settings. 

An initial estimation is made of the relative weighting of 
the two basis attenuation spectra based on a linear combi- 
nation of the basis spectra. Then, a reiterative calculation is 
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performed that converges upon the relative weightings of the 
basis attenuation spectra that produce the observed measure- 
ments. A separate reiterative calculation is performed for 
each voxel. In this way, an accurate measurement of the 
attenuation spectrum for each voxel is provided. 5 

DESCRIPTION OF THE FIGURES 

FIG. 1 (Prior Art) shows a CT imaging device as known 
in the art. 

10 

FIG. 2 shows exemplary attenuation spectra for bone, soft 
tissue and fat. 

FIGS. 3a-3c illustrate an essential restriction on charac- 
teristics of attenuation spectra used in the present method. 

15 

DETAILED DESCRIPTION 

FIG. 1 shows a CT imaging device as known in the art of 
CT imaging. The device includes an X-ray source 20, such 
as an electron tube, and X-ray detectors 22. An object 24 to 
be imaged is disposed between the source 20 and detectors 20 
22. The source and detectors are attached to a rigid frame 26 
that is free to rotate around the object 24. X-ray beam paths 
28 extend from the source to detectors through the object 24. 

In operation, X-ray intensity is measured at every detector at 
many different orientations of the rigid frame with respect to 25 
the object. The number of discrete measurement made in a 
complete scan is equal to NxK, where N is the number of 
positions of the frame used during a scan, and K is the 
number of detectors. In the following discussion, it is 
understood that there are NxK beam paths, with each beam 30 
path corresponding to an X-ray intensity measurement. 

The object is divided into volume elements, or voxels. 
Each voxel within the object is located in the paths of several 
beam paths. Therefore, each voxel affects many X-ray 35 
intensity measurements. The challenge of CT imaging is to 
separate the X-ray attenuation contributions from each 
voxel. 

THEORY 

40 

The attenuation along any beam path is a function of the 
attenuation coefficient of the material in the beam path and 
the energy spectrum of the X-ray source. For a given beam 
path r through the object, the corresponding detector mea- 
surement C is given by: 45 

C(r)=-logJ S(E)exp(J r /i(x,E)dx)dE 

where: 

1) S(E) is the energy spectrum of the X-ray source and must ^ 

be known to the experimenter, and 

2) /i(x,E) is the linear attenuation coefficient of a voxel at 

X-ray photon energy E and location x. 

Consider the operator B defined by B u=C. The operator B 
is nonlinear with respect to fa ^ 

This nonlinear characteristic is the main cause of beam 
hardening artifacts in CT images. 

The problem of calculating CT images while minimizing 60 
beam hardening errors can be restated as follows: given S(E) 
and C(r), find /^(x,E) that satisfies C(r)=-log/S(E)exp(J„w 
(x,E)dx)dE. 

The problem of constructing a CT image without beam 
hardening errors is to find jW(x,E) from the detector mea- 65 
surements C(r). Since /*(x,E) is a three-dimensional 
function, and C(r) is a two dimensional function, the prob- 
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lem is not well posed. In order to find a stable solution, 
^(x,E) must be regularized. The present invention includes 
two different regularization methods associated with two 
different CT methods: single energy reconstruction and dual 
energy reconstruction. 

FIG. 2 shows typical relative attenuation coefficients fa(E) 
for bone, soft tissue and fat. The attenuation coefficients are 
monotonically decreasing with increasing E. 

SINGLE ENERGY RECONSTRUCTION 

In the single energy reconstruction embodiment of the 
present invention, it is necessarily assumed that each voxel 
consists of at most two different known materials. The 
materials may be different in each voxel. For a particular 
voxel at location x, the attenuation coefficient //(x,E) can be 
written as a linear combination of two attenuation coeffi- 
cients according to the equation 

Ax,E)=v(x> 1 (x,E)+(1-v(x)> 2 (x,E) 

where J w 1 (x,E) is the attenuation coefficient for material 1, 
fa(x, E) is the attenuation coefficient for material 2, and v(x) 
is the volume fraction of material 1 within the voxel at 
location x. Other voxels within the object may contain 
materials other than materials 1 and 2, but each voxel is 
modeled to have exactly two distinct materials. A CT image 
is provided by finding the volume fraction v for each voxel. 

In the present invention it is necessary to impose the 
following requirements on the attenuation coefficients fa 
and fa: 

1) Both fa and fa must be monotonically decreasing with 
increasing E, and 

2) The quantity fa-fa must be either positive and mono- 
tonically decreasing, or negative and monotonically 
increasing. 

Although not always specifically noted in the present 
description, it is understood that the attenuation coefficients 
fA are functions of location x and photon energy E, and the 
volume fraction v is a function of location x. 

FIGS. 3a-3c illustrate requirement (2). FIG. 3a shows 
attenuation coefficients which cross, and therefore violate 
requirement (2). FIG. 3b shows attenuation coefficients for 
which fa- fa is either positive mono tonic increasing or 
negative monotonic decreasing and therefore violate 
requirement (2). FIG. 3c shows an exemplary plot of attenu- 
ation coefficients that satisfy requirement (2). 

Next we define an operator P that is linear in /a: 

P iM (r)=J J S(E)/<x,E)dEdx. 

It is noted that P -1 exists and that P -1 P^I. This is true 
because 

J S(E)/z(x,E)dE=P -1 P/z. 

According to the single energy reconstruction method of 
the present invention, the volume fraction v(x) in the voxels 
comprising the object is found using the following algo- 
rithm: 

1) Assume an initial estimate for v(x), given by v°(x). 

2) Calculate a total estimated attenuation coefficient using 
the equation and the current estimate of v(x), given by 
v*(x) (v°(x)=v*(x) in the first iteration): 

/=v*(xK(x, E)+(l-v*(x)K(x,E) 

3) Calculate ju k+1 using the equation: 

/ +1 = p- 1 {P/-B/+C(r)} 
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where C(r) is the measured data. 

4) Solve for v k+1 in the following equation: 

/+ 1 =V*+l^ 1 +(l— V k+1 )fl2 

where: 

// 1 =JS(E) i w 1 (x,E)dE, and 
^=JS(EK(x,E)dE. 

5) Substitute v* +1 for v k in step (2), and reiterate the process 
until successive values v k+2 , v* +3 , v k+4 . . . converge. In 
one particular embodiment, convergence is defined by: 


II V +1 - C(r)\\ 2 
II Bi* - C(r)|| 2 


> 0.98. 


D(x) are found. The attenuation spectrum may be used to 
identify materials within the voxel, or for other purposes. 

According to the dual energy reconstruction method of 
the present invention, the basis coefficients A and D are 
5 found by using the following algorithm: 

1) Assume initial estimates for basis coefficients A(x) and 
D(x), given by A°(x) and D°(x). 

2) Calculate /u k using the equation and the current esti- 
mates of A and D, given by A*(x) and D*(x) (A°(x)=A*(x) 
10 and D°(x)=D*(x) in the first iteration): 

//(x,E)=A a (x)// 1 (x,E)+D a (x) j w 2 (x ; E) 

3) Using /u k , calculate ju k+1 ’ sl for the spectrum S 1 (E) using 
the equation: 

/ + i’ 5 i=p-i{P/-B//+ C ( r )}, 


Generally, the algorithm should be iterated at least twice 
for significant improvement over conventional techniques. 
More preferably, the algorithm is repeated 4-5 times. Sig- 20 
nificant improvements are generally not provided by repeat- 
ing the algorithm more than 5 times. 

The initial estimate v° does not need to be close to the 
final value of v° in order to arrive at a correct solution using 
the present method. However, an initial estimate close to the 25 
final value of the convergent series will reduce the number 
of iterations and the processing time. 

It is noted that the present invention cannot be applied to 
a single voxel in isolation, unless the object comprises a 
single voxel. This is because calculations for every voxel 30 
affect calculations for other voxels, and the entire solution 
must be self-consistent. 

DUAL ENERGY RECONSTRUCTION 

35 

In the dual energy reconstruction, it is assumed that the 
attenuation coefficient //(x,E) can be expressed as a linear 
combination of two basis attenuation functions ^(XjE) and 
ft(x, E): 

4D 

Ax,E)=A(x) i w 1 (x,E)+D(x) i M2(x,E), 

where basis coefficients A(x) and D(x) are real numbers. The 
basis functions i w 1 (x,E) and fa(x, E) do not necessarily cor- 
respond to attenuation coefficients of any known material. In 
a particular object comprising many voxels, each voxel may 45 
be modeled by different basis attenuation functions. For 
example, a single object may comprise several basis attenu- 
ation functions, but any particular voxel is associated with 
exactly two basis attenuation functions. 

Just as in the single energy method, it is necessary to 50 
impose the following requirements on the basis attenuation 
functions and fa\ 

1) Both fa and fa must be monotonically decreasing with 
increasing E, and 

2) The quantity fa- fa must be either positive and mono- 55 
tonically decreasing, or negative and monotonically 
increasing. 

In the dual energy reconstruction, two scans of the object 
must be performed with distinct X-ray beam spectra, S 1 (E) 
and S 2 (E). The beam spectra S 1 (E) and S 2 (E) can be 60 
provided by applying different accelerating voltages to an 
X-ray tube (e.g. 80 kV and 120 kV as known in the art), or 
by placing different X-ray filters in the X-ray beam. Many 
different techniques for altering the X-ray spectrum are well 
known in the art and are often used in X-ray tomography. 65 
In the dual energy reconstruction, the attenuation spec- 
trum for a voxel is characterized when coefficients A(x) and 


where operators P, B, and C(r) are defined with respect to 
S/E). 

4) Using ju k , calculate /u k+1 ’ S2 for the spectrum S 2 (E) using 
the equation: 

/ + i’ 52 =P — i{p/_ B / +C ( r )}, 

where operators P, B, and C(r) are defined with respect to 
S 2 (E). 

5) Calculate A k+1 and D* +1 using the following simultaneous 
equations: 

J^ +1,S1 = A k+1 fa sl +D k+1 fa sl , and 
^^=A^fa S2+ D^fa S2 , 

where // 151 =jS 1 (E) i a 1 (E)dE, i a 152 =jS 2 (E)// 1 (E)dE, fa sl = 
jS 1 (E)// 2 (E)dE, and 

/^2,s2 = jS 2 (E) i w 2 (E)dE . 

6) Substitute A* +1 and D* +1 in step (2), and reiterate the 
process until successive values A k+1 , A k+2 . . . and D* +1 , 
T> k+2 . . . converge. In one particular embodiment, conver- 
gence is defined by: 


II y +1 - C(r)\\ 2 
\W - C(r)|| 2 


> 0.98 for S\(E) and S 2 {E) 


After being computed, the coefficients A and D and basis 
attenuation functions fa and fa can be used to reconstruct the 
attenuation spectra of individual voxels. The reconstructed 
attenuation spectrum can be used to identify materials in the 
object, or for many other purposes such as bone density 
measurements. 

As in the single energy method, the dual energy method 
cannot be applied to a single voxel in isolation, unless the 
object comprises a single voxel. This is because calculations 
for every voxel affects calculations for other voxels, and the 
entire solution must be self-consistent. 

The methods of the present invention provide improve- 
ments in beam hardening errors. Further applications of the 
present invention to signal processing from X-ray tomog- 
raphy will be apparent to one skilled in the art. 

PROOF 

This section contains the theoretical foundation for the 
iterative Single energy reconstruction algorithm. We show 
that the reconstruction problem can be formulated as a fixed 
point problem and the Banach fixed point theorem ensures 
that the algorithm will always converge to a unique solution. 
For ease of understanding, we divide the proof into two 
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theorems. Theorem 1 uses the results from Theorem 2 to 
illustrate the unique convergence properties of the iterative 
algorithm. 

Theorem 1 (Iterative polychromatic reconstruction) 

Let m*(x,E) be an unknown attenuation function which 5 
can be expressed by a two mixture model: 

where v*(x) is the volume fraction of the first material at 
location x and fa, fa are attenuation spectra of the two 
known materials, with the following properties: 

1) For any x, fa and fa are monotonically decreasing 
functions with respect to E. 

2) For any x, the quantity fa-fa is a positive monotonically 

decreasing function with respect to E. 15 

Let C=B/i* be the polychromatic projection of m* with 

respect to the polychromatic spectrum S(E). Then the 
sequence {v*} which satisfies: 

v A ^ 1 +(l-v /: )/i 2 =p- 1 {P/- 1 -B/- 1 +C} 20 

where 

A=JS(EK(x,E)dE, 25 

and 

Ar/S(EK(x,E)dE, 

30 

converges to the true volume fraction function v* indepen- 
dent of the initial estimate v°. 

In the above problem, the true volume fraction v* is a 
fixed point mapping T, 

V*=T(v*), 35 

where T is defined as 


8 

|P 0 (r)-P t (r)|S|E{f}-E{g}HP A ,-P,4 
Taking p -1 on both sides of the above equation leads to 

I J r( w a(x)-w 6 (x))(« 1 (x)-/i 2 (x))dx |<|JXv fl (x)-v 6 (x))C« 1 (x)-/i 2 (x))dx| 

The above inequality is equivalent to 

||w fl — W fe p||v fl — vj 
IT(vJ-T(v 6 )p||v fl -vJ 

where the norm is defined as 

ve V, | |v| |= J | J r v(r) (m, (x)-/i 2 (x))dx|dr 

Thus T is a contraction on V. 

Theorem 2 

Let f and g be positive monotonically decreasing func- 
tions such that f — g is a positive monotonically decreasing 
function, then 

E{f}-E{g}>[-log E {exp(-f)}]-[-log E{exp(-g)}] 

where the expectation is taken with respect to any probabil- 
ity density function. 

Proof 

Let the expectation be based on probability density func- 
tion h. We define t:[0,l]-^9Jas 

t(a)=E A {g+a(f-g)}-[-log E A {exp(-g-a(f-g))}] 

By Jensens inequality, t(0)^0. Differentiating with 
respect to a, we have 

J - g)exp(-g - g(f - g))} 

da~ h 8 E h {exp(-g-a(f -g))} 

^=E h \f-g}-E h ,{f-g} 


w =T(v): w// 1 +(l-w> 2 =P 1 {P(v^i+(1-v> 2 )-B(v^ 1 +(1-v)// 2 )+C} 

The Banach fixed point theorem states that if T is a 40 
contraction on the space of volume fraction V, 

Vx,ye V,\ae5I,a ^ 1 such that | |T(x)— T(y) 1 1 = ct| |x— y | ^ 

then T has precisely one fixed point and the iterative 
sequence {v*}, from the procedure v*=T(v* -1 ) converges to 45 
the unique fixed point v* of T with arbitrary v°eV. Hence, 
proving that T is a contraction on V is sufficient. 

Let v fl , v b be any two volume fraction functions and fa, fa, 
p a , p b , w fl , w b be defined as follows: 

50 


where h' is a probability density function defined by 

h , = exp(-g - a(f - g)) 

E h {exp(-g - a(f - g))} 

Since f and g are monotonically decreasing functions, 
exp( — g — a(f — g — )) is a monotonically increasing func- 
tion for any ae[0,l]. 

Thus we conclude that 


— = E h {f -g}- E h ,{f - g} > 0 for any a e [0, 1] 


Mb=VtJh+Q--Vb)M2 

p a =PjLl a -BjLl a +C 

p b =P/* b -B/t b +C 

w fl =T(v fl ) 

w 6 =T(v 6 ) 

Defining f(r,E)= J^(x,E)dx and g(r,E)= J/^(x,E)dx, we 
have 

|P a (r)-P 6 (r)|=|E{f }— E{g}— [— log E{exp(-f)}]+[-log E{exp(-g)}]| 

where the expectation is taken with respect to the beam 
spectrum profile S(E). Using the results from the above 
equation leads to 


This shows that t is a monotonically increasing function. 
Thus we have t(l)>t(0), which is equivalent to: 

E{f}-E{g}>[-l°g E{exp(-f)}]-[-log E{exp(-g)}]. 

It will be clear to one skilled in the art that the above 
embodiment may be altered in many ways without departing 
from the scope of the invention. Accordingly, the scope of 
60 the invention should be determined by the following claims 
and their legal equivalents. 

What is claimed is: 

1. A method for calculating characteristics of an object 
comprising voxels using an X-ray computer assisted tomog- 
65 raphy device, comprising the steps of: 

a) projecting through the object an X-ray beam having a 
known energy spectrum indicated by S(E); 
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b) detecting the X-rays with an X-ray detector; 

c) repeating steps (a) and (b) for at least two different 
orientations of the X-ray beam with respect to the 
object; 

d) defining a first material and a second material within 5 
each voxel, wherein the first and second materials have 
known distinct attenuation spectra, given by J w 1 (x,E) 
and /^(x, E), respectively, such that a total attenuation 
coefficient //(x,E) of the object is given by: 

//(x,E)=v(x> 1 (x,E)+(1-v(x))/< 2 (x,E) 


10 

g) defining a first basis attenuation function and a second 
basis attenuation function, given by i a 1 (x,E) and ju 2 (x, 
E), respectively, such that a total attenuation coefficient 
/u(x, E) of the voxel is given by: 

/f(x,E)= A(x )// 1 (x, E) +D (x)// 2 (x, E) 

where coefficients A and D are real numbers; 

h) defining initial coefficients A* and D*; 

i) calculating an estimated total attenuation coefficient 
//(x,E) from the equation; 


where v(x) is a number between 0 and 1 and indicates a 
volume fraction of each voxel comprised of the first mate- 
rial; 15 

e) defining an initial estimated volume fraction v*(x); 

f) calculating an estimated total attenuation coefficient 
ju\x,E) from the equation; 

/(x ; E)=v a (xK(x,E) + (1-v*(x)K(x,E) 20 


g) calculating a new estimated total attenuation coefficient 
/u k+1 (x,E) from the following equation: 


/ + 1 (x,E)=p- 1 {P/(x,E)-B/(x,E)+C(r)} 25 

where P is an operator defined by V i u=\ r \ S(E)/i(E)dEdx, B 
is an operator defined by B//=-logjS(E)exp(J^w(E)dx)dE, 
and C(r) represents X-ray detector output, and r indicates the 
X-ray beam path; 

h) calculating a new estimated volume fraction v* +1 from 30 
the following equation: 


/ + 1 (x,E)=v Xc+ 1 (x,E)/i 1 (x,E)+(l-v ife+ 1 (x))/i 2 (x,E) 

where /^(E)* /S(E) i w 1 (E)dE and i w 2 (E)=jS(E) i w 2 (E)dE 35 

i) repeating steps (f), (g) and (h) using the new estimated 
volume fraction v* +1 (x) instead of the initial estimated 
volume fraction v*(x). 

2. The method of claim 1 wherein steps (f), (g), and (h) are 
repeated using a successive value of the estimated volume 40 
fraction v k+2 . 

3. The method of claim 1 wherein steps (f), (g), and (h) are 
repeated 4 times using successive values of the estimated 
volume fraction. 

4. The method of claim 1 wherein steps (f), (g), and (h) are 45 
repeated until 


11 V +1 - cwiiz 

11 ^ - C(r)\\ 2 


> 0.98. 


50 


5. A method for calculating characteristics of an object 
comprising voxels using an X-ray computer assisted tomog- 
raphy device, comprising the steps of: 

a) projecting through the object a first X-ray beam having 

a known energy spectrum indicated by S 1 (E); 55 

b) detecting the first X-ray beam with an X-ray detector; 

c) repeating steps (a) and (b) for at least two different 
orientations of the first X-ray beam with respect to the 

object; 60 

d) projecting through the object a second X-ray beam 
having a known energy spectrum indicated by S 2 (E); 

e) detecting the second X-ray beam with an X-ray detec- 
tor; 

f) repeating steps (d) and (e) for at least two different 65 
orientations of the second X-ray beam with respect to 
the object; 


/(x)(x,E)=A a (x)// 1 (x,E)+D a (x)// 2 (x,E) 


j) calculating a new estimated total attenuation coefficient 
^r+i,5i( x ,g) associated with spectrum S 1 (E) from the 
following equation: 

// + 1 ’‘ sl (x,E)=P 51 _ 1 {P 51 j « A (x,E)-B 51 j « A (x,E)+C 1 (r)} 

where p« is an operator defined by P 53 //=J r jS 1 (E) i w(E) 
dEdx, B S1 is an operator defined by B 51i a=-logjS :L (E)exp 
(J r (E)dx)dE, and C-^r) represents X-ray detector output 
associated with S 1 (E); 

k) calculating a new estimated total attenuation coefficient 
// + ’ 52 (x,E) from the following equation: 

// + 1 ’ 52 (x,E)=P 52 - 1 {P 52 « A (x ; E)-B 52 « A (x ; E) + C 52 (r)} 

where P52 is an operator defined by P 52 a=J r /S 2 (E) i w(E) 
dEdx, B 52 is an operator defined by B 52 a=-log jS 2 (E)exp 
( J„w(E)dx)dE, and C 2 (r) represents X-ray detector output 
associated with S 2 (E); 

l) calculate new estimated coefficients A* +1 (x) and D* +1 
(x) from the following simultaneous equations: 

// +1 ’‘ s1 (x,E)=A a+1 (x)// 1iS1 (x,E)+D a+1 (x) j w 2 51 (x,E), and 

/ +1 ’ 52 (x,E)=A^ 1 (x> 1 ^ 2 (x,E) + D" +1 (xK^ 2 (x,E), 


where 


kji( x )= J S 1 (EK(x,E)dE, k,s2( x )=J S 2 (EK(x,E)dE, 
Jksi ( X W S 1 (E) J M 2 (x,E)dE, and M2,s 2 ( x )=I S 2 (E)/i 2 (x,E)dE; 


m) repeating steps (i), (j), (k) and (1) using the new 
estimated coefficients A* +1 (x) and D* +1 (x) instead of 
the initial estimated coefficients A*(x) and D*(x). 

6. The method of claim 5 wherein steps (i), (j), (k) and (1) 
are repeated using successive values of the estimated coef- 
ficients A k+2 (x) and D* +2 (x). 

7. The method of claim 5 wherein steps (i), (j), (k) and (1) 
are repeated 4 times using successive values of the estimated 
coefficients. 

8. The method of claim 5 wherein steps (i), (j), (k) and (1) 
are repeated until 


\\Bsifi N -c^h 


> 0.98 and 


lifted 1 - C 2 (r)|| 2 
II W* - c 2 (r)|| 2 


£ 0.98, 


where N is an integer indicating number of algorithm 
iterations. 





